Gravitational waves from pulsations of neutron stars 
described by realistic Equations of State 
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In this work we discuss the time-evolution of nonspherical perturbations of a nonrotating neutron 
star described by a realistic Equation of State (EOS). We analyze 10 different EOS for a large sample 
of neutron star models. Various kind of generic initial data are evolved and the gravitational signals 
are computed. We focus on the dynamical excitation of fluid and spacetime modes and extract the 
corresponding frequencies. We employ a constrained numerical algorithm based on standard finite 
differencing schemes which permits stable and long term evolutions. Our code provides accurate 
waveforms and allows to capture, via Fourier analysis, the frequencies of the fluid modes with an 
accuracy comparable to that of frequency domain calculations. The results we present here are 
useful for providing comparisons with simulations of nonlinear oscillations of (rotating) neutron star 
models as well as testbeds for 3D nonlinear codes. 
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I. INTRODUCTION 

Neutron stars (NSs) are very compact stars that are 
born as the result of gravitational collapse [1] . They are 
highly relativistic objects and their internal composition, 
governed by strong interactions, is, at present, largely 
unknown. After NS formation (either as the product of 
gravitational collapse or of the merger of a binary NS 
system), nonisotropic oscillations are typically present. 
These oscillations are damped because of the emission of 
gravitational waves (GWs). In general, the nonspheri- 
cal oscillations of a NS are characterized by two types of 
proper modes (called quasi-normal modes, QNMs here- 
after): fluid modes, which have a Newtonian counterpart, 
and spacetime (or curvature) modes, which exist only in 
relativistic stars and are weakly coupled to matter. See 
Ref. [2, 3] for a review. The QNMs frequencies carry 
information about the internal composition of the star 
and, once detected, they could, in principle, be used to 
put constraints on the values of mass and radius and thus 
on the Equation of State (EOS) of a NS [4]. 

In principle, only 3D simulations in full nonlinear Gen- 
eral Relativity (GR) with the inclusion of realistic mod- 
els for the matter composition (as well as electromagnetic 
fields) can properly investigate the neutron star birth and 
evolution scenarios. The first successful steps in this di- 
rection have been recently done by different groups [5— 
8] . However the complexity of the physical details behind 
the system and the huge technical/computational costs of 
these simulations are still not completely accessible, and 
alternative/approximate approaches to the problem are 
still meaningful. In particular, we recall the work of Dim- 
melmeier et al. [9] , who simulated oscillating and rotating 
NSs (described by a polytropic EOS) in the conformally 
flat (CF) approximation to GR and using specific initial 



data. Another approximate, and historically important, 
route to studying NS oscillations is given by perturbation 
theory, i.e. by linearizing Einstein's equation around a 
fixed background (see Ref. [10]). The perturbative ap- 
proach has proved to be a very reliable method to under- 
stand the oscillatory properties of NS as well as a useful 
tool to calibrate GR nonlinear numerical codes. 

Although most of the work in perturbation theory has 
been done (and is still done) using a frequency-domain 
approach (in order to accurately compute mode frequen- 
cies), time-domain simulations are also needed to com- 
pute full waveforms [11-23]. In particular, Allen et 
al. [11], via a multipolar expansion, derived the equa- 
tions for the even parity perturbations of spherically sym- 
metric relativistic stars and produced explicit waveforms. 
They argue that, for various kinds of initial data, both 
fluid and spacetime modes are present, but they have dif- 
ferent relative amplitudes depending on the initial excita- 
tion of the system. Technically, the problem is reduced to 
the solution of a set of 3 wave- like hyperbolic equations, 
coupled to the Hamiltonian constraint: two equations for 
the metric variables, in the interior and in the exterior of 
the star, and one equation for the fluid variable in the in- 
terior. The Hamiltonian constraint is preserved (modulo 
numerical errors) during the evolution. Ruoff [12] derived 
the same set of equations directly from the ADM [24] 
formulation of Einstein's equations and used a similar 
procedure for their solution. This work showed that the 
presence of spacetime modes in the waveforms strongly 
depends on the initial data used to initialize the evo- 
lution. In particular, conformally flat initial data can 
totally suppress the presence of spacetime modes. Both 
studies use a simplified description of the internal com- 
position of the star, i.e. a polytropic EOS with adiabatic 
exponent V — 2. In addition Ref. [12] explored also the 
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use of one realistic EOS. The author found a numerical 
instability related to the dip in the sound speed at neu- 
tron drip point. This instability was independent either 
of the formulation of the equations or of the numerical 
(finite differencing) scheme used. The use of a particular 
radial coordinate was proposed to cure the problem. 

In this work we reexamine the problem of the evolu- 
tion of the perturbation equations for relativistic stars 
investigating systematically the gravitational radiation 
emitted from the oscillations of nonrotating neutron star 
models described by a large sample of realistic EOS. We 
use, specified to the Regge- Wheeler gauge and a static 
background, the general gauge-invariant and coordinate- 
independent formalism developed in [25-28] . The result- 
ing system of equations is equivalent to the formulation 
of [11, 12]. For the even-parity perturbation equations, 
we adopt a constrained numerical scheme [13, 15, 23], 
(different from any of those adopted in [11, 12]) which 
permits long-term, accurate and stable evolutions. We 
use standard Schwarzschild-like coordinate system and 
we don't need the technical complications of Ref. [12]. 
We evolve various kind of initial data (for odd and even- 
parity perturbations) for 47 neutron star models com- 
puted from 10 different EOS. We compute and show grav- 
itational waveforms and extract QNMs frequencies. Our 
accurate results show that there are no relevant qualita- 
tive differences in the waveforms with respect to previous 
work limited to polytropic EOS. This was expected since, 
in first approximation, the features of the waves depend 
only on the star mass and radius (in particular on the 
compactness). The results we report are comprehensive 
data obtained with a new, efficient numerical code and 
they complete the information already present in the lit- 
erature. 

The plan of the paper is as follows. In Sec. II we briefly 
review the formalism used and the equations describing 
the nonspherical perturbations of a spherically symmet- 
ric star. In Sec. Ill and Sec. IV the construction of the 
equilibrium star models and the EOS sample are dis- 
cussed. Sec. V deals with the initial data setup, and 
in Sec. VI we present the results. We use dimensionless 
units c = G = Mq = 1, unless otherwise specified for 
clarity purposes. 



try, are decoupled 1 . 

In this work we assume the Regge- Wheeler gauge [29] . 
The background metric of a static spherical star 
of radius R and mass M, obtained by solving the 
Tolman-Oppenheimer-Volkoff (TOV) equations [24] of 
hydrostatic equilibrium (sec below), is written in 
Schwarzschild-like coordinates as 

g^dx^dx" = -e 2a dt 2 +e 2/3 dr 2 +r 2 (d8 2 +sm 2 9d(t) 2 ), (1) 

where a(r) and f3(r) are function of r only. The matter 
is modeled by a perfect fluid: 



= {e+p)u»u v + pg< MJ , 



(2) 



where p is the pressure, u* 1 the fluid 4-velocity, and e = 
p(l + e) the total energy density. Here p denotes the rest- 
mass density and e the specific internal energy. The rest- 
mass density can also be written in terms of barionic mass 
tub and the baryonic number density n as p = msn. The 
speed of sound is defined as C 2 = dp/de. The adiabatic 
exponent is 
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In the Regge- Wheeler gauge, the even-parity metric per- 
turbation multipoles are parametrized by three (gauge- 
invariant) scalar functions (k, \, ip) as 2 
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where Yi m are the usual scalar spherical harmonics. Here, 
k is the perturbed conformal factor, while x is the actual 
GW degree of freedom. Since the background is static, 
the third function tp is not independent from the others, 
but can be obtained from k and x solving the equation 
(for r<R) [27] 
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II. PERTURBATION EQUATIONS 

The perturbation equations are obtained by specializ- 
ing to the nonrotating case the general gauge-invariant 
and coordinate-independent formalism for metric pertur- 
bations of spherically symmetric spacetimes introduced 
by Gerlach and Sengupta [25, 26] and further developed 
by Gundlach and Martin-Garcia [27, 28]. Let us recall 
that, due to the isotropy of the background spacctime, 
the metric perturbations can be decomposed in multi- 
poles, i.e. expanded in tensorial spherical harmonics. 
These are divided in axial (or odd-parity) and polar (or 
even-parity) modes which, due to the spherical symmc- 



where the mass function m = m(r) is defined as e~ 2/3 ( r ) = 
1 — 2m(r)/r and represents the mass of the star inside a 
sphere of radius r. This equation also holds for r > R 
with m(r) = M. 

In addition, when the background is static the metric 
perturbations are actually described by two degrees of 
freedom, (k,x), only in the interior [27], while only one 
degree of freedom remains in the exterior. A priori there 



1 Under a parity transformation ((9, <f>) —> (tt — 8, tt + </>)) the axial 
modes transform as (— and the polar modes as (— 1)*. 
We omit hereafter the multipolar indexes for convenience of no- 
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is no unique way of selecting which evolution equations 
to use for numerical simulations (the ones most conve- 
nient mathematically could not be so numerically) so that 
different formulations of the problem have been numeri- 
cally explored in the literature [11-13, 27]. In particular, 
Rcf. [13] showed that it can be useful to formulate the 
even-parity perturbations problem using a constrained 
scheme, with one elliptic and two hyperbolic (wave-like) 
equations. One hyperbolic equation is used to evolve % 
in the interior and exterior; the other hyperbolic equa- 
tion serves to evolve, in the interior, the perturbation of 
the relativistic enthalpy H = Sp/(p + e), where 5p is the 
pressure perturbation. The system is closed by the ellip- 
tic equation, the Hamiltonian constraint, that is solved 
for k. Following Ref. [13] we express the equations in 
term of an auxiliary variable S = x/r, whose amplitude 
tends to a constant for r — > oo and thus is more conve- 
nient for the numerical implementation. We recall that 
the variable S is the same used by Ruoff [12] and the 
relationship with the variables of Allen ct al. [11] is given 
by k = Fx\\cn/ r an d S = e~ 2a SAiien- In the star interior, 
r < R, the evolution equation for S reads 



S. tt + e 2 ^S, rr = e 2a \- 



+ 



4nr(5p — e) + 



6m 



r 2 \ r ) r 2 



/ TfX \ ~ i ' / 1 1 

2e {^2 +4:7Tr P) + 8n£ - — 



6m 



(6) 



the one for H becomes 

-H tt + Cy^H rr = e 2 » { [^(1 + Cl) (7) 
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where A = £(i + 1). Eqs. (6) and (8) are also valid in 
the exterior, with C 2 = p = H = 0, m(r) = M and 



e 2a = 1 — 2M/r. Since in the star exterior the spacetime 
is described by the Schwarzschild metric, the perturba- 
tion equations can be combined together in the Zerilli 
equation [30] 



*!« + ^ (c) * (c) = o 



(9) 



for a single, gauge- invariant, master function \I>( e ), the 
Zerilli-Moncrief function [30, 31]. The function is 
the Zerilli potential (see for example [32]) and r* = 
r + 2M ln[r/ (2M) - 1] is the Regge- Wheeler tortoise co- 
ordinate. In terms of the gauge-invariant functions \ an d 
k, reads 



\]/( e ) 



2r(r - 2M) 
A[(A-2)r + 6M] 
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(10) 



The inverse equations can be found, for instance, in 
Ref. [12]. In our notation they read 
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Let's mention briefly the boundary conditions to impose 
to these equations. At the center of the star all the func- 
tion must be regular, and this leads to the conditions: 



X ~ r 
k 
H 



£+2 



(13) 
(14) 
(15) 



At the star surface S is continuous as well as its first 
and second radial derivatives. On the contrary, k and its 
first radial derivative are continuos but k^ rr can have a 
discontinuity due the term 8n(p + e)H/C 2 in Eq. (8). At 
the star surface, r = R, Eq. (7) reduces to an ODE for 
H, that is solved accordingly. 

On a static background, the odd-parity perturbations 
are described by a single, gauge-invariant, dynamical 
variable that is totally decoupled from matter. This 
function satisfies a wave- like equation of the form [25, 33] 



(o) 



with a potential 
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This equation has been conveniently written in terms of 
the "star-tortoise" coordinate r* defined as df^/dr = 



4 



e 13 ^ 01 . In the exterior, f* reduces to the Regge- Wheeler 
tortoise coordinate r* introduced above and Eq. (16) be- 
comes the well-known Regge- Wheeler equation [29] . The 
relation between \I/( ) and the odd-parity metric multi- 
poles is given, for example, by Eqs. (19)-(20) of Ref. [32]. 

The principal quantities we want to obtain are the 
gauge-invariant functions vf( e /°). These functions are di- 
rectly related to the "plus" and "cross" polarization am- 
plitudes of the GWs by (see e.g. [32, 34]): 



h + ~ ih x = \ jr £ N t (*W + i$ £>) _ 2 Y tm (e, 0), 

' =2 m=-e 



where N e = y/(£ + 2)(t + l)t(£ - 1) and 
spin- weighted spherical harmonics of spin- weight s 
The GWs luminosity at infinity is given by 



(18) 

vYim are the 
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TABLE I: A list of the EOS name and references that we use 
in this work. 



Name 


Authors 


References 


A 


PflnnnflTiTiflflnp 


[OOJ 


B 


Pandharipande 


f36l 


c 


Bethe and Johnson 




FPS 


Lorenz, Ravenhall and Pethick 


[38, 39] 


G 


Canuto and Chitre 


[40] 


L 


Pandharipande and Smith 


[41] 


N 


Walecka and Serot 


[42] 


O 


Bowers, Gleeson and Pedigo 


[43, 44] 


SLy 


Douchin and Haensel 


[45] 


WFF 


Wiringa, Fiks and Farbroncini 


[46] 



-=-£y>, 



£=2 m= 



+ 



(19) 



where the overdot stands for derivative with respect to 
coordinate time t. The energy spectrum reads 



dE 

dbJ 167T 2 



V 



£=2 m= 



(20) 



where indicates the Fourier transform of V P^°\ 

uo = 2-kv and v is the frequency. 



III. EQUILIBRIUM STELLAR MODELS 

The equilibrium configuration of a spherically symmet- 
ric and relativistic star is the solution of the TOV equa- 
tions 



m, r = 4irr e , 
P,r = -{s + p)a ir , 
1 (m + Anr 3 p) 
a ' r ~ 2 (r 2 - 2mr) 

with the boundary conditions 

m(0) = 
p(R) = 



(21) 



a(R) 



In 1 



2M 



(22) 
(23) 

(24) 



The system is closed with an EOS p(p). Eq. (23) formally 
defines the star radius, R. Eq. (21) define the structure 
of the fluid and its spacetime in the interior, i.e. r < R, 
then the solution is matched at the exterior Schwarzschild 
solution, Eq. (24). 



IV. EQUATIONS OF STATE 

Neutron stars are composed by high density baryonic 
matter. The exact nature of the internal structure, deter- 
mined essential by strong interactions, is unknown. To 
model the neutron star interior (approximated) many- 
body theories with effective Hamiltonians are usually em- 
ployed. The principal assumptions are that the matter is 
strongly degenerate and that it is at the thermodynam- 
ics equilibrium. Consequently, temperature effects can 
be neglected and the matter is in its ground state (cold 
catalyzed matter). Under this conditions the EOS has 
one-parameter character: e(n) and p(n) or e(p). 

The composition of a neutron star consists qualita- 
tively of three parts separated by transition points, see 
Fig. 1. At densities below the neutron drip, e < e d ~ 10 11 
gem -3 (the outer crust) the nuclei are immersed in an 
electron gas and the electron pressure is the principal 
contribute to the EOS. In the inner crust, £d < e < 10 14 
gem -3 , the gas is also composed of a fraction of neu- 
trons unbounded from the nuclei and the EOS softens 
due to the attractive long-range behaviour of the strong 
interactions. For e > 10 14 gcm~ 3 a homogeneous plasma 
of nuclcons, electrons, muons and other baryonic matter 
(e.g. hyperons), composes the core of the star. In this re- 
gion the EOS stiffens because of the repulsive short-range 
character of the strong interactions. The bottom panel 
of Fig. 1 shows the adiabatic exponent which, under the 
assumption of thermodynamics equilibrium, determines 
the response of pressure to a local perturbation of den- 
sity. We mention that, as explained in detail in Ref [45], 
in a pulsating NS the "actual adiabatic exponent" can 
be higher than that obtained from the EOS because the 
timescale of the beta processes are longer than the dy- 
namical timescales of the pulsations. Thus the fraction 
of particles in a perturbed fluid clement is assumed fixed 
to the unperturbed values (frozen composition). We re- 
fer the reader to [47] for all the details on neutron star 
structure and the complex physics behind it (e.g. elastic- 
ity of the crust, possible supcrfluid interior core, magnetic 
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FIG. 1: Pressure (top) and the adiabatic exponent (bottom) 
as a function of the total energy density for various EOS. 
Notice here we are using cgs units. 



fields). 

Most of the mass (see the bottom panel of Fig. 2) is 
constituted by high density matter, so that the maximum 
mass of the star is essentially determined by the EOS of 
the core. On the other hand, the star radius strongly 
depends on the properties of the matter at low densities, 
due to Eq. (23). 

To compensate the ignorance on the interior part of the 
star it is common to consider a large set of EOS derived 
from different models. In our work we employ 7 realistic 
EOS already used in [41], and in many other works (see 
for example Refs. [4, 10, 48-50] on pulsations of relativis- 
tic stars, and [51-53] on equilibrium models of rotating 
stars). Maintaining the same notation of [41], they are 
called A, B, C, G, L, N and O EOS. Most of the models 
in the sample are based on non-relativistic interactions 




5 5.5 6 6.5 7 7.5 8 8.5 9 9.5 
R 




FIG. 2: Top panel: mass versus radius for the NS models of 
Table II. Bottom panel: profile of the total energy density e 
for all models with M = 1.4. 



modeled with Reid soft core type potentials. EOS N [42] 
and O [43, 44] are instead based on relativistic interaction 
and many-body theories. Model G [40] is an extremely 
soft EOS, while L [41] is extremely stiff. EOS A [35] and 
C [37] are of intermediate stiffness. In addition, we use 
the FPS EOS [38, 39], and the SLy EOS [45], modeled by 
Skyrme effective interactions. The FPS EOS, in particu- 
lar is a modern version of Friedman and Pandharipandc 
EOS [39]. The last EOS considered is the UV14+TNI 
(here renamed WFF) EOS of [46], which is an intermedi- 
ate stiffness EOS based on two-body Urbana UV14 po- 
tential with the phenomenological three-nucleon TNI in- 
teraction. The composition is assumed to be of neutrons. 
For all the EOS models the inner crust is described by 
the BBS [54] or the HP94 [55] EOS, while for outer crust 
the BPS EOS [56] is used. We refer to Table I and cited 
references for further details. 

Realistic EOS are usually given through tables. To use 
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TABLE II 


: Neutron Star models. 


From left to right the 


columns report: 


the name 


of the model, the EOS type, 


the mass M, 


the radius 


R, the compactness M/R, the central energy 


density e c and the central 


pressure p c . 




Model 


EOS 


M 


R 


M/R 


£c 


Pc 


A10 


A 


1.00 


6.55 


0.15 


1.96xl0 -3 


2.35xl0 -4 


A12 


A 


1.20 


6.51 


0.18 


2.38xl0~ 3 


3.76xl0~ 4 


A14 


A 


1.40 


6.39 


0.22 


3.01 xl0~ 3 


6.50xl0~ 4 


A16 


A 


1.60 


6.04 


0.26 


4.46xl0~ 3 


1.49xl0" 3 


Amx 


A 


1.65 


5.60 


0.29 


6.78xl0~ 3 


3.24xl0~ 3 


BIO 


B 


1.00 


5.78 


0.17 


1.37xl0~ 3 


5.13xl0~ 4 


B12 


B 


1.20 


5.53 


0.22 


1.59xl0~ 3 


1.02xl0" 3 


B14 


B 


1.40 


4.96 


0.28 


1.88xl0~ 3 


3.42xl0~ 3 


Bmx 


B 


1.41 


4.73 


0.30 


4.65xl0~ 3 


5.33xl0~ 3 


CIO 


C 


1.00 


8.15 


0.12 


3.35xl0~ 3 


1.05xl0~ 4 


C12 


c 


1.20 


8.05 


0.15 


4.45xl0~ 3 


1.67xl0~ 4 


C14 


c 


1.40 


7.90 


0.18 


7.76xl0~ 3 


2.65xl0~ 4 


C16 


c 


1.60 


7.67 


0.21 


9.76xl0~ 3 


4.49xl0~ 4 


Cmx 


c 


1.85 


6.69 


0.28 


1.18xl0~ 3 


1.94xl0" 3 


FPS10 


FPS 


1.00 


7.32 


0.14 


1.44xl0~ 3 


1.50xl0~ 4 


FPS12 


FPS 


1.20 


7.30 


0.16 


1.75xl0~ 3 


2.30xl0~ 4 


FPS14 


FPS 


1.40 


7.22 


0.19 


2.22xl0~ 3 


3.64xl0~ 4 


FPS16 


FPS 


1.60 


7.02 


0.23 


4.77xl0~ 3 


6.38xl0~ 4 


FPSmx 


FPS 


1.80 


6.22 


0.29 


1.45xl0~ 3 


2.50xl0~ 3 


G10 


G 


1.01 


5.76 


0.17 


1.73xl0~ 3 


5.36xl0~ 4 


G12 


G 


1.20 


5.42 


0.22 


2.10xl0~ 3 


1.20xl0~ 3 


Gmx 


G 


1.36 


4.62 


0.29 


2.72xl0~ 3 


5.65xl0~ 3 


L10 


L 


1.00 


9.61 


0.10 


5.55xl0~ 3 


4.15xl0~ 5 


L12 


L 


1.20 


9.76 


0.12 


3.52xl0~ 3 


5.60 xl0~ 5 


L14 


L 


1.40 


9.89 


0.14 


5.03xl0~ 3 


7.41xl0~ 5 


L16 


L 


1.60 


9.98 


0.16 


7.81 xl0~ 4 


9.75xl0~ 5 


Lmx 


L 


2.68 


9.23 


0.29 


5.80xl0~ 4 


9.17xl0~ 4 


N10 


N 


1.00 


8.81 


0.11 


6.38xl0~ 4 


5.70xl0~ 5 


N12 


N 


1.20 


8.98 


0.13 


7.05xl0~ 4 


7.59xl0~ 5 


N14 


N 


1.40 


9.13 


0.15 


7.81 xl0~ 4 


9.97xl0~ 5 


N16 


N 


1.60 


9.24 


0.17 


2.33xl0~ 3 


1.30xl0~ 4 


Nmx 


N 


2.63 


8.65 


0.30 


7.08 xl0~ 4 


1.20xl0" 3 


O10 





1.00 


8.29 


0.12 


1.41xl0~ 3 


7.39xl0~ 5 


012 





1.20 


8.41 


0.14 


1.64xl0~ 3 


1.03xl0~ 4 


014 





1.40 


8.50 


0.16 


1.96xl0~ 3 


1.43xl0~ 4 


016 





1.60 


8.54 


0.19 


2.45xl0~ 3 


1.97xl0~ 4 


Omx 


o 


2.38 


7.75 


0.31 


5.15xl0~ 3 


1.66xl0" 3 


SLylO 


SLy4 


1.00 


7.78 


0.13 


7.72xl0~ 4 


1.12xl0~ 4 


SLyl2 


SLy4 


1.20 


7.80 


0.15 


8.41 xl0~ 4 


1.66xl0~ 4 


SLyl4 


SLy4 


1.40 


7.78 


0.18 


9.17xl0~ 4 


2.45xl0~ 4 


SLyl6 


SLy4 


1.60 


7.70 


0.21 


2.57xl0~ 3 


3.70xl0~ 4 


SLymx 


SLy4 


2.05 


6.71 


0.30 


8.74xl0~ 4 


2.50xl0~ 3 


WFF10 


WFF3 


1.00 


7.29 


0.14 


9.86xl0~ 4 


1.46xl0~ 4 


WFF12 


WFF3 


1.20 


7.30 


0.16 


l.lOxlO -3 


2.18xl0~ 4 


WFF14 


WFF3 


1.40 


7.26 


0.19 


1.23xl0~ 3 


3.34xl0~ 4 


WFF16 


WFF3 


1.60 


7.14 


0.22 


3.36xl0~ 3 


5.48xl0~ 4 


WFFmx 


WFF3 


1.84 


6.40 


0.29 


1.18xl0~ 3 


2.18xl0~ 3 
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them in a numerical context it is necessary to interpolate 
between the tabulated values. The interpolation can not 
be chosen arbitrarily but must properly take into account 
the First Law of Thermodynamics, see e.g. [57], that, in 
the case of a temperature independent EOS, reads: 

A thermodynamically consistent procedure is described 
in [58] and it is based on Hermite poynomials. Essentially 
the method permits to interpolate a function forcing the 
match on the tabulated points both of the function and 
of its derivatives. We implement this scheme using cubic 
Hermite poynomials as already done in [51], The proce- 
dure used is described in detail in Appendix A. 

V. INITIAL DATA 

In principle the choice of the initial data for the per- 
turbation equations should take into account, at least 
approximately, the astrophysical scenario in which the 
neutron star is born. Such scenario could be, for example 
the gravitational collapse or the merger of two neutron 
stars. Only long-term simulations in full general relativ- 
ity can investigate highly nonlinear and nonisotropic sys- 
tem until they settle down in a, almost spherical, quasi- 
cquilibrium configuration. A pcrturbative analysis, like 
the one we propose, could then start from this point once 
the metric and matter tensors had been projected along 
the corresponding (tensorial) spherical harmonics. This 
approach, in principle possible, is however beyond the 
scopes of the present work. 

Inspired by previous perturbative calculations [11, 12], 
we consider different kinds of initial data such that they 
are the simplest, well-posed and involve perturbations of 
both the fluid and/or the metric quantities. 

In the case of the even parity perturbations we start 
the evolutions from 3 different initial excitations of fluid 
and matter variables: 

1. Conformally Flat Initial Data. We set 5(0, r) = 
and give a fluid perturbation of type: 

H(0,r)=A^y~\m(ir(n + l)^). (26) 

The function k(0,r) is computed consistently solv- 
ing the Hamiltonian constraint. The profile of H in 
Eq. (26) is chosen in order to approximate the be- 
havior of an enthalpy eigenfunction with n nodes. 
In this way, only some modes can be (prominently) 
excited. Since we would like to focus on the princi- 
pal fluid modes we chose a zero-nodes initial data 
setting n = 0. 

2. Radiative Initial Data. We set fc(0,r) = and 
H(0,r) as Eq. (26). The function 5(0, r) is com- 
puted consistently solving the Hamiltonian con- 
straint. 



3. Scattering-like Initial Data. We set H (0, r) = 
and *(°)(0,r) as a Gaussian pulse: 

^ c )(0,r)^Aexp(^- (r ~ 2 rc)2 ^, (27) 

with r c = 70 M and b = M. The functions 
fc(0, r) and 5(0, r) are computed consistently from 
Eqs. (11)-(12). 

Initial data of type 1 and 2 are chooscn to be time- 
symmetric (S t t — k,t = H t = 0). On the one hand, 
this choice can be physically questionable because the 
system has an unspecified amount of incoming radiation 
in the past. On the other hand, it is the simplest choice 
and guaratees that the momentum constraints are triv- 
ially satisfied and only the Hamiltonian constraint needs 
to be used for the setup. The Gaussian in type 3 initial 
data is ingoing, i.e. \& t = $ rt , and the derivatives of 
the other variables are computed consistently with this 
choice. 

In the case of odd-parity perturbations we start the 
evolution using an ingoing narrow gaussian as in Eq. (27) 
for *(°)(0, r). The amplitude of the perturbation A is 
everywhere chosen equal to 0.01. 

VI. RESULTS 

For each EOS, we study a rapresentative set of mod- 
els with M = 1, 1.2, 1.4, 1.6 and a model whose mass 
is close to the maximum mass allowed, for a total of 
47 neutron stars. The principal equilibrium properties 
are summarized in Table II. Fig. 2 shows the mass- 
radius diagram for all the models computed. The star 
radius spans a range from R ~ 5 (in the case of EOS 
B and G) to R ~ 9 - 10 (for EOS L and N). The or- 
der of stiffness of the EOS can be estimated, on average, 
as: G<B<A<FPS<WFF<SLy<C<0<N<L. For all the 
models described by a particular EOS, the compactness 
M/R increases from about 0.1 to 0.3, corresponding to 
the increase of the star mass and the decrease of the ra- 
dius. The (total) energy density profile as function of the 
radial coordinate in the bottom panel of Fig. 2 explains, 
as discussed above, that most of the mass is due to mat- 
ter with density comparable to the central (maximum) 
density of the star. 

For each model, we numerically evolve the equations 
for the odd and even-parity perturbations described in 
Sec. II using the initial data presented in Sec. V. All the 
details of the numerical schemes employed can be found 
in Appendix B. In the following sections we will discuss 
the results obtained focusing on the I = 2 (quadrupolc) 
multipole, as this is the principal responsible of the grav- 
itational wave emission. The waves are extracted at dif- 
ferent radii, r obs = [50, 100, 200, 300]M. We checked the 
convergence of the waves and the differencies between the 
extraction at r obs = 200M and r obs = 300M are very 
small, so that we can infer to be sufficiently far away 
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from the source. The gravitational waveforms we discuss A. Axial Waveforms 

in the following have always been extracted at the far- 
thest observer, r obs = 300M and they are plotted versus 

observer's retarded time u = t — r° hs . The gravitational waveform that results from scatter- 

ing of Gaussian pulses of GWs off the odd-parity poten- 
tial exhibits the well known structure (precursor-burst- 
ring down-tail [59]) analogue to the black holes case (see 
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for example Ref. [60] for the case of polytropic EOS). 
The characteristic signature of the star in the waveform 
is contained in the ringdown part, which is shaped by 
high frequencies, (quickly) exponentially damped oscilla- 
tions: the w- modes [61]. These modes are pure spacetime 
vibrations and are the analogue of black hole QNMs for 
relativistic stars [2, 62]. 



As a representative case, because the global qualitative 
features are common to all EOS, Fig. 3 exhibits wave- 
forms computed only with EOS A. The compactness of 
the model increases from top to bottom; the left panels 
exhibit the waveforms on a linear scale, while the right 
panels their absolute values on a logarithmic scale. Since 
the damping time increases with the star compactness, 
the maximum mass model, Amx, presents the longest w- 
mode ringdown. On the contrary, model A10 exhibits 
only a one-cycle, small-amplitude ringdown oscillation 
that quickly disappears in the power-law tail. 



In principle, an analysis of the frequency content of the 
axial waveforms by looking at the Fourier spectra or by 
means of a fit procedure based on a quasi-normal mode 
template is possible. Let us focus on model Amx, that 
presents the longest and clearest ringdown waveform. Us- 
ing the same fit-analysis method of [60] we estimate the 

frequency of the fundamental w-modes to be v$ = 9452 
Hz, with a damping time r4°^ ~ 0.07 ms. For comparison, 
we note that a Schwarzschild black hole of the same mass 
has the fundamental I = 2 frequency and damping time 
equal to, respectively, v Bli = 7317 Hz and r BH = 0.09 
ms. We have also computed the energy spectrum of the 
waveform starting from u ~ 170 (the first zero after the 
burst): the spectrum has a single peak centered at a fre- 
quency that differs from uffl of about a few percents. 
However, as discussed in [60], we found that this infor- 
mation is in general very difficult to extract, especially 
for the lowest mass models, due to the rapid damping of 
the modes and their localization in a narrow time win- 
dow. The comparison with frequency domain data (see 
Ref. [60] for polytropic EOS and the discussion in Ap- 
pendix B for realistic EOS) shows that the errors on the 
numbers presented above are of the order of 5%. The 
error on the frequencies increases up to about 12% for 
models with M ~ 1.4 and to about 20% for models with 
smaller mass. Moreover, the damping times can not be 
reliably estimated with a fit procedure when they are 
too short (see Table III in Appendix B). In summary, 
although the analysis of the waveforms through a fit pro- 
cedure works for some particular models, in general it 
seems uncapable to give numbers as robust and reliable 
as those provided by a standard frequency domain ap- 
proach. See for example Ref. [49] for details about this 
approach. 



x lCf 3 
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100 200 300 400 500 
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FIG. 4: Even parity waveforms for model WFF14 generated 
by initial data of type 1 (solid line) and type 2 (dashed line). 
In the first case, only fluid modes (the /-mode and the first 
p-mode) are present; in the second case, 10-mode oscillations 
are present at early times (0 < u < 50). The inset shows 
long-term evolution (corresponding to a total time of ~ 20 
ms) used to compute the energy spectrum of Fig. 5. 



B. Polar Waveforms 

We start this section presenting together waveforms 
generated by conformally flat (type 1) and radiative (or 
non-conformally flat, type 2) initial data. Former studies 
with polytropic EOS [12, 23] showed that initial data of 
type 1 determine the excitation of fluid modes only. On 
the contrary, initial data of type 2 produce a gravitational 
wave signal where both spacetime and fluid modes are 
present. 

As expected, this qualitative picture is confirmed also 
for realistic EOS. Figure 4 exhibits the waveform \I>( e ) for 
the representative model WFF14. For conformally flat 
initial data (solid line) the figure shows that the Zcrilli- 
Moncrief function oscillates at (mainly) one frequency, of 
the order of the kHz. The corresponding energy spectrum 
(solid line in Fig. 5) reveals that the signal is in fact 
dominated by the frequency Vf — 2126 Hz, but there is 
also a second peak at v p = 6909 Hz. This two frequencies 
are recognized as those of the fundamental fluid mode / 
and of the first (pressure) p-mode. 

For non-conformally flat initial data (dashed line in 
Fig. 4) the first part of the signal, i.e. < u < 50, is 
dominated by a high-frequency and strongly-damped os- 
cillation typical of curvature modes. For u > 50, the type 
1 and type 2 waveforms are practically superposed. The 
corresponding energy spectrum (dashed line in Fig. 5 has, 
superposed to the two narrow peaks of the fluid modes, a 
wide peak centered at higher frequency (~ 11 kHz) that 
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V [Hz] v [Hz] 



FIG. 5: Energy spectra of the type 1 (solid line) and type 2 
(dashed line) initial data evolution of Fig. 4. The two narrow 
peaks at 2126 Hz and 6909 Hz correspond to the /-mode and 
the first p-mode frequencies. The wide peak at ~ 11 kHz 
corresponds to u>-mode excitation. 



is typical of the presence of spacetime excitation [2] . 3 

The information of Fig. 5 is complemented by the 
Fourier spectrum of the metric variable S in Fig. 6. In 
contrast to the fluid variable H , which contains only nar- 
row peaks for both kind of initial data, for type 2 wave- 
forms the metric variable S also exhibits a broad peak 
which is absent for type 1 initial data. Note, that the 
picture that we have discussed so far for model WFF14 
remains qualitatively unchanged for all the other EOS. 

Since our numerical scheme allows us to evolve the 
system in time as long as we wish, we can produce very 
long time-series to accurately extract, via Fourier analy- 
sis, the fluid mode frequencies. We did this analysis sys- 
tematically for all the models considered. In Tables IV- 
V in Appendix B we list the frequencies of the /-mode 
and the first p-mode extracted from the energy spec- 
tra. By comparison with the published frequencies of 
Andersson and Kokkotas [4] for some models with EOS 
A (obtained via frequency domain calculations), we esti- 
mate that the errors on our values are tipically smaller 
than 1%. For a fixed EOS, the frequencies increase with 
the star compactness. For a model of given mass, the 
/-mode frequency generally decreases if the EOS stiff- 
ens. The same (on average) is true for the first p-mode 



Note that, in order to obtain the cleanest fluid-mode peaks, in 
the Fourier transform of type 1 waveform we discarded the first 
four GW cycles, which are contaminated by a transient due to 
the initial excitation of the system. On the other hand, for type 2 
waveform, we considered the full time-series. In this case, it is not 
possible to cut the first part of the signal because it also contains 
the curvature mode contribution that we want to analyze. 



FIG. 6: Power spectrum of the variable S for model WFF14 
in the case of initial data of type 1 and 2. 



frequency. Following Ref. [4], we present in Fig. 7 the 
frequencies that we ha ve computed as a function of the 
mean density y/M/R 3 (/-mode) and of the compactness 
M/R (p-mode) of the star. Globally, they show a very 
good quantitative agreement with previously published 
results calculated by means of a standard frequency do- 
main approach [4]. 

We conclude this section by discussing the wave- 
forms generated by initial data of type 3. This kind of 
"scattering-type" initial condition constitutes the even- 
parity analogue of that discussed in Sec. VI A for the odd- 
parity case. In Fig. 8 we show the waveforms from 
3 models of EOS WFF. The waveforms in the top panel 
of the figure are very similar to those of the left panels of 
Fig. 3. The logarithmic scale (bottom panel) highlights 
the main qualitative difference: i.e., fluid- mode oscilla- 
tions are present, in place of the nonoscillatory tail, after 
the w-mode ringdown. Note that, in principle, the tail 
will emerge in the signal after that all the fluid modes 
have damped (i.e., on a time scale of a few seconds). 

The process of w-mode excitation is instead exactly 
the same as for the odd parity case: the ringdown phase 
is longer (and thus clearly visible) for the more compact 
models. The frequencies are also very similar. For exam- 
ple, for model WFFmx (the one discussed in the figure) 

we have v w — 8638 Hz and a damping time Tw ~ 0.05 
ms, while for model Amx vffl = 9798 Hz and rffl c± 0.05 
ms (to be compared with v$ = 9452 and t4°^ — 0.07 
ms). 



VII. CONCLUSIONS 

In this work we have discussed the time-evolution of 
nonsphcrical (matter and gravitational) perturbations of 
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FIG. 7: Comparison between the oscillation properties of 
models computed with different EOS. Top panel: /-mode fre- 
quencies as a function of y?M/R 3 . Bott om panel: the first 
p-mode frequency (multiplied by the mass M) as a function 
of the compactness M/R. 



FIG. 8: Excitation of even-parity w-modes for some WFF 
models from given initial data of type 3. As in the odd-parity 
case (compare with Fig. 3), the ringdown is more pronounced 
for more compact models. 



nonrotating neutron stars described by a large sample 
of realistic EOS. The current study extends the work of 
Allen et al. [11] and Ruoff [12], who focused essentially 
on polytropic, but one, EOS models. We have used an 
improved version of a recently developed ID perturbative 
code [13] that has been thoroughly tested and used in the 
literature [15, 16, 23]. 

The main, new result, presented here is that our con- 
strained numerical scheme allows us to stably evolve the 
even-parity perturbation equations without introducing 
any "special" coordinate change, as it was necessary in 
Ref. [12]. In addition, despite the EOS that we con- 
sider are very different, the outcome of our computa- 
tions is fully consistent (as expected) with previous stud- 
ies involving polytropic EOS. In particular: (i) for even- 
parity perturbations, if the initial configuration involves 



a fluid excitation, (type 1 and 2 initial data), the Zerilli- 
Moncrief function presents oscillations at about 2-3 kHz 
due to the excitation of the fluid QNMs of the star; 
(ii) if we set S ^ at t = (type 2, i.e. the non- 
conformally flat condition is imposed), high frequencies, 
strongly damped w-mode oscillations are always present 
in the waveforms; (iii) the w-mode excitation is gener- 
ally weak, but it is less weak the more compact the star 
model is; consistently, (iv) for scattering-like initial data 
in both the odd and even-parity case the presence of w- 
modes is more striking the higher is the compactness of 
the star 4 ; even-parity fluid modes are typically weakly 



4 Note that type 3 initial data are non conformally flat as well, 
since \ ^ 0. 
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excited in this case. 



APPENDIX A: TABLE INTERPOLATION 



Thanks to the long-term and accurate evolutions we 
can perform, we extracted the fluid mode frequencies 
from the Fourier transform of time series of the waves 
with an accuracy comparable to that of frequency domain 
codes. For what concern the frequencies and damping 
times of spacetime modes, pretty good estimates can be 
obtained when damping times are not too short; i.e., for 
the more compact models. When these frequencies will 
be revealed in a gravitational wave signal, they will hope- 
fully provide useful information on the internal structure 
of neutron stars. In particular, recognizing both fluid 
and w-modes in the signal could permit, in principle, to 
estimate the values of mass and radius of the NS and thus 
to put strong constraints on the EOS model [4, 48, 63]. 
We have limited our analysis to the first two fluid modes 
because they are the most responsible for gravitational 
waves emission. We checked that, by changing the ini- 
tial fluid perturbation, our evolutionary description per- 
mits also to easily capture the frequencies of higher over- 
tones [15]. 

In addition, despite all the approximation that we have 
introduced (initial data, no rotation, no magnetic fields), 
we believe that the approach to NS oscillations described 
in this paper can provide physical information comple- 
mentary to that available from GR nonlinear evolutions. 
In particular, it can also be used to provide useful test- 
beds for GR nonlinear codes and must be seen as a first 
step to compare/constrast with nonlinear simulations of 
neutron star oscillations. In conclusion, since we can 
freely specify the initial data of the metric and matter 
variables at initial time, one can also think to use our 
present tool to evolve further in time an almost spherical 
configuration that is the outcome of a long-term numer- 
ical relativity simulation. A perturbative evolution like 
the one we discussed here (possibly complemented by a 
complex 3D (magneto)-hydrodynamics source, as an im- 
provement of the approach discussed in Ref. [13]) could 
then start when the 3D fully nonlinear simulation ends. 
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Here we describe the method used to interpolate the 
tables of the EOS. The interpolation scheme is based on 
Hermite polynomials and was introduced in Ref. [58] . We 
follow [51] by using 3rd order (cubic) polynomials and we 
list here the relevant formulas for completeness. Consider 
a function y(x), given the tables of yj and y'j and a point 
xj < x < Xj+i the interpolated value y(x) is: 



y(x) = y l H {w) + y i+1 H (l - w) 



(Al) 



where 



Axj = {x 



w{x) = 



Axj 



Xj), 



and 



H (w) = 2w 3 - 3w 2 + 1, 
Hi(w) = w 3 - 2w 2 + w 



(A2) 
(A3) 



(A4) 
(A5) 



are the cubic Hermite functions. The principal properties 
of this method arc that for x — > Xs 



2/0) 

y'{x) 



V'r 



(A6) 
(A7) 



The EOS are usually given as 3 column tables with 
the values for n, e and p. In this work we need to com- 
pute e(p) and the speed of sound C 2 . The thermody- 
namics consistency can be achived by computing, for a 
given value of p, first n(p) and then e(n) imposing the 
derivative through Eq. (25). Finally the speed of sound is 
computed consistently with the interpolation. To obtain 
more accurate numerical data, we perform such calcu- 
lations not directly on the functions, but taking loga- 
rithms. Violation of the first law of thermodynamics are 
typically less than 0.1%. We tested also other interpola- 
tion schemes, i.e. linear and spline interpolation, that in 
general gives violation of the thermodynamic principle of 
some percents. 

As a check of the implementation, we evolved a poly- 
tropic model both with the analytic EOS and with ta- 
bles of different numbers of entries. The results perfectly 
agree and we did not find any dependence on the tabu- 
lated points. In addition, the use of other interpolation 
schemes did not produce significant differences: this fact 
suggests that the global numerical errors of the code are 
dominant over the errors related to the violation of the 
thermodynamics principle. In the case of some tables 
(EOS FPS, SLy4 and L), we found that "high-order" in- 
terpolation (cubic Hermite and spline) did not permit an 
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accurate reconstruction of the sound speed. This was 
due to spurious oscillations introduced by the high order 
derivatives. As a consequence, the code gave unphysical 
results. In all these cases, we adopted the linear interpo- 
lation. 



APPENDIX B: NUMERICAL DETAILS, CODE 
TESTS AND MODE FREQUENCIES 

The code we employ in this work is a development 
of that described in [23] and successfully used in many 
works [13, 15-17]. 

The TOV equations (21) are integrated numerically 
(from the center outward), for a given central pressure 
p c (see Table II), using a standard fourth-order Runge- 
Kutta integration scheme with adaptive step size. 

To evolve numerically the perturbations equations, we 
introduce an evenly spaced grid in r with uniform spac- 
ing Ar and we adopt finite differencing approximation 
schemes for the derivatives. In particular, in the con- 
struction of the the computational grid the origin r = 
is excluded and the first point is located at r = Ar/2. 
The resolution is measured as the number of point Ji in- 
side the star radius. The star surface is located at a cell 
center R = Ar/2 + (Ji - l)Ar. 

The hyperbolic evolution equations for S and H are 
all solved with standard, second order convergent in time 
and space, leapfrog algorithm. For the evolution of the 
even-parity equations, the Hamiltonian constraint is used 
to update, at every time step, the variable k. This ellip- 
tic equation is discretized in space at second order and 
reduced to tridiagonal linear system that is inverted. For 
this reason the evolution scheme of the polar perturba- 
tions can be considered a constrained evolution. The 
inner boundary conditions of Eqs. (13)-(15) are imple- 
mented by setting to zero the variables at the first grid 
point. At the outer boundary, standard radiative Som- 
merfcld conditions are imposed. 

The Zerilli-Moncrief function has been obtained in two 
(independent) ways. On the one hand, it has been com- 
puted from S (x) an d k using Eq. (10) for every value of r. 
On the other hand, it has been computed using Eq. (10) 
only at the star surface (r = R, the matching point) and 
then evolved using Eq. (9). The second method has been 
used only as an independent consistency check and all 
the results discussed in this paper are obtained using the 
first method. 

In some situations, we found the convergence of the 
Zerilli-Moncrief computed from Eq. (10) particularly del- 
icate. For example, in case of type 3 initial data (an 
even-parity Gaussian pulse of gravitational radiation) we 
realized that a very accurate computation of the radial 
derivative k. r is needed to obtain an accurate and reli- 
able from k and S. Fig. 9 summarizes the kind of 
problem that one can find computing the Zerilli-Moncrief 
function too naively. It refers to ty( e >(r) at t = 0. We 
consider (for simplicity of discussion) a polytropic model, 
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FIG. 9: Convergence of the Zerilli function for initial data of 
type 3 for a polytropic EOS model. See text for discussion. 



TABLE III: Frequencies v and damping times r of fluid and 
spacetime modes of some models of Andersson and Kokko- 
tas [4], described by EOS A, computed via our approach. 
The frequencies are expressed in Hz and the damping times 
in ms. For the sake of comparison, we report also the values 
of their Table A.l with the "AK" superscript. 



M 




4 K v P 


u AK 


AK 


■ 


1.653 


3080 


3090 7825 


7838 


9342 9824 


0.062 0.064 


1.447 


2580 


2579 7843 


7818 


10038 11444 


0.057 0.027 


1.050 


2183 


2203 7555 


7543 


11267 14328 


0.059 0.017 


i.e. p 


= Kp T 


with K = 


100, r 


= 2 and p c 


1.28 x 10~ 3 



with type 3 initial data (i.e. a Gaussian pulse with a = M 
centered at r = 40). We fix by Eq. (27), we compute 
X and k from Eqs. (11)-(12), we compute k yT . numerically 
and then we we reconstruct ty( e ) via Eq. (10). Fig. 9 
shows that, for low resolution (Ji = 100), using a second- 
order finite-differencing standard stencil to compute fc r is 
clearly not enough, as the reconstructed (thicker dashed- 
dot line) and the "exact" (solid line) are very dif- 
ferent. Incrasing the resolution (to 500 points) improves 
the agreement, which is however not perfect yet (thin- 
ner dashed-dot line). A visible improvement is obtained 
using higher order finite-differencing operators: the fig- 
ure shows that a 4th-order operator (already in the low- 
resolution case with Ji — 100 points) is sufficient to have 
an accurate reconstruction of the Zerilli-Moncrief func- 
tion. The conclusion is that one needs to use at least 
4th-order finite-differencing operators to compute accu- 
rately \I>( e ) from x and k. If this is not done, the resulting 
function is not reliable and it can't be considered a so- 
lution of Eq. (9). Typically, we have seen that, when 
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TABLE IV: Frequencies of the fluid modes for EOS A, B, C, 
FPS and G. From left to right the columns report: the name 
of the model, the EOS type, the /-mode frequency and the 
first p-mode frequency. 



Model 


EOS 


v s [Hz] 


Up [Hz] 


A10 


A 


2146 


7444 


A12 


A 


2305 


7816 


A14 


A 


2512 


7873 


A16 


A 


2833 


7823 


Amx 


A 


3107 


7848 


BIO 


B 


2705 


7440 


B12 


B 


3044 


8061 


B14 


B 


3577 


8893 


Bmx 


B 


3732 


9091 


CIO 


C 


1617 


5052 


C12 


C 


1802 


5349 


C14 


c 


1933 


5606 


C16 


c 


2114 


5876 


Cmx 


c 


2627 


6421 


FPS10 


FPS 


1884 


6326 


FPS12 


FPS 


2028 


6590 


FPS14 


FPS 


2174 


6764 


FPS16 


FPS 


2325 


6891 


FPSmx 


FPS 


2820 


7031 


G10 


G 


2750 


7510 


G12 


G 


3097 


8164 


Gmx 


G 


3891 


9080 



this kind of inaccuracy is present, the amplitude of (part 
of) the Zerilli-Moncrief function (usually the one related 
to the w-modc burst) grows linearly with r instead of 
tending to a constant value for r — > oo. 

We performed extensive simulations to test the code. 
The scheme is stable and permits us to accurately evolve 
the equations as long as we wish. To check convergence of 
the waves, we run different resolutions and we computed 
the energy emitted at infinity (last detector) integrating 
separately the 1 = 2 odd and even contribute of Eq. (19) 
over the evolution time interval. The value of the energy 
converges correctly up to second order terms, 0(Ar 2 ). 

To validate the physical results of the code we com- 
pare the frequencies extracted from our simulations with 
values computed via a frequency domain approach. In 
the case of the polytropic EOS, the results reported in 
Ref. [15] showed that the errors on fluid frequencies are 
less than 1%. Spacetime frequencies for polytropic EOS 
models were checked in Ref. [60]. In this case, when 
the damping times are sufficiently long (i.e. when a nar- 
row Gaussian pulse is used and the star model is very 
compact), it is possibile to estimate v and t using a fit 
procedure, with an error of the order of 6%. 

The accuracy of the frequencies does not change when 
we use realistic EOS. To validate this assertion, we com- 



TABLE V: Frequencies of the fluid modes for EOS L, N, O, 
SLy and WFF. From left to right the columns report: the 
name of the model, the EOS type, the /-mode frequency and 
the first p-mode frequency. 



Model 


EOS 


Vf [Hz] 


Up [Hz] 


L10 


L 


1217 


4599 


L12 


L 


1297 


4850 


L14 


L 


1353 


5025 


L16 


L 


1395 


5115 


Lmx 


L 


1871 


5109 


N10 


N 


1415 


5324 


N12 


N 


1466 


5694 


N14 


N 


1497 


5893 


N16 


N 


1522 


5960 


Nmx 


N 


1952 


5470 


O10 


N 


1481 


5924 


012 


N 


1578 


6200 


014 


N 


1643 


6283 


016 


N 


1734 


6217 


Omx 


N 


2217 


5969 


SLylO 


SLy 


1691 


5818 


SLyl2 


SLy 


1804 


6087 


SLyl4 


SLy 


1932 


6279 


SLyl6 


SLy 


2029 


6468 


SLymx 


SLy 


2607 


6601 


WFF 


WFF 10 


1889 


6544 


WFF 


WFF 12 


1973 


6766 


WFF 


WFF 14 


2126 


6909 


WFF 


WFF 16 


2282 


7016 


WFF 


WFFmx 


2718 


7016 



pared the frequencies extracted from our waveforms with 
those of Andersson and Kokkotas [4] for EOS A with 
mass M = 1.653, M = 1.447 and M = 1.050 (see Ta- 
ble A.l of Ref. [4]). Our results are listed in Table III, 
together with the data of [4] for completeness. Fluid fre- 
quencies are typically captured with an accuracy below 
1%, while spacetime frequencies and damping times can 
be estimate with decent accuracy (5%) only for the max- 
imum mass model. As a consequence, we expect that a 
similar accuracy for fluid modes, i.e. of the order of 1%, 
should be expected for the 47 NS models of Table II. For 
completeness, the corresponding frequencies are listed in 
Table IV and Table V. 

The simulations we have discussed in the main text 
of the paper use resolutions of Jj = 400 and J; = 800 
respectively for the odd and even-parity case. The 
Courant-Friedrichs-Lewy factor is set to At/Ar = 0.4. 
For each model, the outer boundary of the grid is at 
r = 600M. The final evolution time is i cnd = 3000M 
for even-parity evolutions and t end = 800M for the odd- 
parity ones. 
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